% this code takes the solution from Gensys and generates Table 4

clearvars -except wd
cd(fullfile(wd, 'Table4')); 

Table_4 = zeros(20,5);

for i_case = 1:5
    filename = ['temp' num2str(i_case) '.mat'];
                
    load(filename) 
    
    % the next few lines just regenerate position indices and n_k weights that 
    % were not saved
    n_sectors = 341;
    
    params.rho_a = 0.9;         % persistence agg tech shocks
    params.rho_mu = 0.9;         % persistence monetary shock
    params.rho_k = 0.9;         % persistence idiosync tech shock
    params.beta = 0.9975;      % discount factor
    params.theta = 6;          % elast of subst within sectors
    params.delta = 0.5;
    params.eta = 2;
    
    pos_ps = zeros(n_sectors,1);
    pos_p = zeros(n_sectors,1);
    pos_Ep = zeros(n_sectors,1);
    pos_mc = zeros(n_sectors,1);
    pos_w = zeros(n_sectors,1);
    pos_yk = zeros(n_sectors,1);
    pos_zk = zeros(n_sectors,1);
    pos_lk = zeros(n_sectors,1);
    pos_ak = zeros(n_sectors,1);
    
    % Position indices for the variables
    pos_c = 1;
    pos_Ec = 2;
    pos_pc = 3;
    pos_Epc = 4;
    pos_mu = 5;
    pos_y = 6;
    pos_a = 7;
    pos_i = 8;
    maxind = 8;
    
    for i = 1:n_sectors
        % superscript k prices
        pos_ps(i) = maxind + i;
        % subscript k variables
        pos_p(i) = maxind + n_sectors + i;
        pos_Ep(i) = maxind + 2*n_sectors + i;
        pos_mc(i) = maxind + 3*n_sectors + i;
        pos_w(i) = maxind + 4*n_sectors + i;
        pos_yk(i) = maxind + 5*n_sectors + i;
        pos_zk(i) = maxind + 6*n_sectors + i;
        pos_lk(i) = maxind + 7*n_sectors + i;
        pos_ak(i) = maxind + 8*n_sectors + i;
    end
    
    Ck = csvread("2002CkDetailed.csv");
    Znum = csvread("2002RevShareDetailed.csv");
    n = length(Ck); % Number of Industries
    Z = Znum(1:n,1:n);
    
    Wc = Ck/sum(Ck); % Share of consumption; the weights omega
    W = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega
    psi = params.delta*(params.theta - 1)/params.theta;
    psi_C = (1-params.delta*(params.theta - 1)/params.theta);
    n_k = inv(eye(n)-psi.*W)*(psi_C*Wc);
    
    % Y = C + Z
    % note Zk is the total intermediate demand by sector k
    % recover C time series
    Cnew = (n_k'*yt_imp_mu(pos_yk,:)  - sum(psi*(n_k'.*W)*yt_imp_mu(pos_zk,:)) )/psi_C
    
    for k = 1:n_sectors
        Ckt(k,:) = yt_imp_mu(pos_c,:) - params.eta*(yt_imp_mu(pos_p(k),:)-yt_imp_mu(pos_pc,:));
    end
    
    % cumulate over time each sector's consumption response
    cum_Ck = sum(Ckt')';
    
    [ZZZ, i_zzz] = sort(cum_Ck);
    Table_4(1:10,i_case) = ZZZ([end-9:end]);
    Table_4(11:20,i_case) = ZZZ([1:10]);
end

format bank
% print out Table 4 results
Table_4

cd(fullfile(wd, 'Table4')); 

% save as files
save('Table4.mat')
filename = ['Table4.csv'];
dlmwrite(filename, Table_4, 'precision', '%.10f', 'delimiter', '\t');
